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We apply the microscopic coupled-cluster method (CCM) 
to the spin-i XXZ models on both the one-dimensional chain 
and the two-dimensional square lattice. Based on a system- 
atic approximation scheme of the CCM developed by us pre- 
viously, we carry out high-order ab initio calculations using 
computer-algebraic techniques. The ground-state properties 
of the models are obtained with high accuracy as functions 
of the anisotropy parameter. Furthermore, our CCM analy- 
sis enables us to study their quantum critical behavior in a 
systematic and unbiased manner. 

PACS numbers: 75.10.Jm, 75.30.Gw, 75.50.Ee 



Over the last few years the microscopic coupled-cluster 
method (CCM) 0] has successfully been applied to many 
lattice Hamiltonian systems ||-|^ , producing the best or 
among the best results in terms of accuracy and power. 
For example, in the so-called SUB2 approximation dis- 
cussed below, the exact analytic solution of the corre- 
sponding CCM equations often exhibits a terminating 
point. This is clearly demonstrated |^ to correspond 
to a physical critical point of the system by the behavior 
within the same approximation of such calculated quanti- 
ties as the order parameter, correlation functions, and the 
gap in the excitation spectrum. Recently we have devel- 
oped several efficient, systematic approximation schemes 
of the CCM specifically tailored for use with lattice sys- 
terns 

In this article, we report our new results from an ah 
initio high-order calculation for the spin-i XXZ model 
on both the one-dimensional chain (IDC) and the two- 
dimensional square lattice (2DSL). In particular, we have 
obtained with high accuracy the ground-state energy, 
the anisotropic susceptibility (i.e., the second-order en- 
ergy derivative with respect to the anisotropy parame- 
ter), and the staggered magnetization, as functions of the 
anisotropy parameter. Furthermore, since these physical 
quantities are obtained as functions of the anisotropy pa- 
rameter, we are able to study the possible quantum phase 
transitions of the anisotropic models in a systematic, un- 
biased manner. This is in contrast to the series expansion 
technique Q in which one has to make a Fade approxi- 
mation or to assume a particular critical behaviour (from 
the spin- wave theory of Anderson 0, for example) for 
those physical quantities from the outset. 



Since the details of our CCM analysis have been pub- 
lished elsewhere we only outline the specific approx- 
imation schemes employed here. The spin system un- 
der consideration is the so-called spin-i XXZ model de- 
scribed by the following Hamiltonian 



H = 



1=1 p=i 



A.sfsf+^ + -(s+s,+^ + .Sj 



(1) 



where the index I runs over all N {-^ oo) lattice sites 
with the usual periodic boundary condition imposed; the 
index p runs over all z nearest-neighbour sites; the oper- 
ators sf and sf (= sf ±?sf) are spin operators; and A is 
the anisotropy parameter. The special case A = 1 gives 
the isotropic Heisenberg model, which, with spin s = ^ 
on the 2DSL, has been under intensive study over the 
last six years or so Q . 

We first consider the case of A — > oo. Equation (1) 
then reduces to the Ising model with a classical ground 
state (i.e., the Neel state) given by two alternating sub- 
lattices, one with all spins down, the other with all spins 
up. For clarity, we use the index {?} exclusively for the 
spin-down sublattice, and the index {j} exclusively for 
the spin-up sublattice. Naturally, we choose the Neel 
state as the model state 1$) in our CCM analysis, and in- 
corporate the quantum correlation effects by considering 
the excitations with respect to this model state. The el- 
ementary operators of these excitations are clearly given 
by the spin-raising operators on the i-sublattice and 
the spin-lowering operators s'J on the j-sublattice. The 
CCM ansatz for the ground ket-state is therefore given 
by 
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with the correlation operators 5*271 defined by 
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where we have restricted ourselves to the conserved sector 
of zero z-component of total spin, s^^^^j (= ■^f )' by 

including only those configurations with equal numbers 
of spin-flips on both sublattices. 
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The ground-state energy and the c-number coefficients 
{Sii...i„,ji...j„} of Eq. (3) are determined by taking the 
inner products of the Schrodinger equation in the form 
e~^He^\^) = Eg\^), firstly with the model state itself 
and secondly with the states constructed by acting on 
1$) with the corresponding correlation operators in S'2„. 
We thus find respectively the ground-state energy Eg, 

Eg = me-'^He^l^); (4) 

and the coupled set of equations for {Si-^...i„,ji...j„}, 

i^K^^. ■ ■ ■ ■ ■ ■ <e-^i/e^|<i>) = 0, (5) 

withn= l,2,...,7V/2. 

Each of the above equations always involves the Hamil- 
tonian in a similarity-transformed form, namely, 

e"^i/e^ = H + [H,S] + ^[[H, S],S] + ---, (6) 

where the expansion series terminates at the fourth order 
1^. Therefore, once an approximation for S is chosen, 
no further approximation is necessary in order to obtain 
and solve the coupled set of equations (5). For the spin-i 
XX Z model of Eq. (1), it is easy to derive the following 
exact equation for the ground state energy per spin, 

f = -|(A + 2M, (7) 

where bi = 5i^,;+p is the nearest-neighbour pair corre- 
lation coefficient, and z — 2,4 for the IDC and 2DSL 
models respectively. We note that bi is independent of 
both the index i and index p by the lattice symmetries. 

We clearly need an approximation method to truncate 
S for any practical calculation. The three most com- 
monly used truncation methods are: the SUBn scheme 
in which all correlations involving only n or fewer spins 
are retained; the simpler SUBn-m sub-approximation 
scheme where only SUBn correlations spanning a range 
of no more than m adjacent lattice sites are included; 
and finally the systematic local LSUBm scheme, which 
includes all possible multi-spin correlations over a spec- 
ified locale on the lattice, where m is the nominal index 
that characterizes the size of the given locale. In each 
case the remaining correlation coefficients are set to zero. 
For example, the LSUB4 scheme for the ID spin-i model 
retains three independent configurations, represented by 
61,63, and 34 respectively In particular, bi corre- 
sponds to the nearest-neighbour two-spin-flip configura- 
tion mentioned before, 63 to the third-nearest-neighbour 
two-spin-flip configuration, and 34 to the spin-flip config- 
uration of four adjacent spins. The corresponding LSUB4 
coupled equations are given by j^], 

1 - 2A6i - 36? + 26163 + 2bl + 2.94 = 0, (8) 
6? - 4A63 - 46163 + 34 = 0, (9) 
-A(6? + 26163) + (74(A + 46i + 63) + 2616^ = 0. (10) 



After solving these coupled equations, we obtain the 
ground-state energy by substituting 61 into Eq. (7). 

For the higher-order approximations the derivation of 
the coupled equations becomes very tedious. We have 
developed our own software using C-l — h and Fortran to 
automate this process. Also, we have used standard 
computer algebra packages to check our results indepen- 
dently. For the LSUBm schemes, we have derived and 
solved the coupled equations up to m = 10 for the IDC 
case and up to m = 6 for the 2DSL case. It should be 
noted that all of our calculations were done on micro- 
computers. The numbers of independent spin-flip config- 
uration coefficients retained in the these two cases are 81 
and 72, respectively. 

Some of our results for the IDC model have already 
been published In particular, we have shown that for 
a given value of m the LSUBm scheme reproduces ex- 
actly the corresponding 2mth order of large-A perturba- 
tion theory, and that the LSUBm scheme also gives good 
results in the planar region (|A| < 1) where perturbation 
theory is not valid. The new high-order calculations fur- 
ther push the numerical results closer to their exact coun- 
terparts obtained by the Bethe ansatz over a wide 
range of A (0 < A < 00). In particular, at the isotropic 
point (A = 1) the LSUBIO scheme yields -0.4420 for 
the ground-state energy per particle. We have made a 
naive attempt to extrapolate our results for the LSUBm 
scheme with m = 4,6,8,10, and find that a 1/m^ rule 
seems to fit them very well. The ground-state energy per 
particle after this extrapolation yields —0.4431 ± 0.0001 
from a least squares fit, while the exact result by Bethe 
ansatz ^ is —0.4432 to the accuracy of four significant 
figures. 

In Fig. 1 we show some of our results for the ground- 
state energies of the 2DSL model, together with a Monte 
Carlo calculation [|l^ on an 8 x 8 lattice for compari- 
son. One sees that the results for our high-order calcu- 
lations are in excellent agreement with the Monte Carlo 
calculations over a wide range. At the isotropic point 
(A = 1), our best numerical results for the ground-state 
energy of the 2DSL model, obtained from the LSUB6 
scheme, is —0.6670. We find that the 1/m^ rule also fits 
our 2D LSUBm data well. The extrapolated result is 
—0.6691 ± 0.0003. This number compares well with the 
values —0.66934 ±0.00004 from a large-scale Monte Carlo 
calculation by Runge 0, and -0.6694±0.0001 from the 
series expansion techniques |6| . 

The more interesting feature from our 2DSL LSUBm 
calculations is the appearance of terminating points in 
the real solutions of the LSUBm equations for m > 2, as 
indicated in Fig. 1, beyond which the solution becomes 
complex. This behavior is totally different from that 
observed in the IDC model where, unlike in the SUB2 
scheme discussed above which retains two-spin correla- 
tions of arbitrarily long range, there is no evidence of 
such terminating points in the localized LSUBm scheme 
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at any value m < 10. Since we have derived the corre- 
sponding coupled equations for a given LSUBm scheme 
as closed analytical forms with A as a parameter, we can 
straightforwardly study the terminating points by taking 
the derivatives on both sides of the coupled set of equa- 
tions (e.g., Eqs. (8-10)) with respect to A and solving 
directly for the derivatives of the coefficients. From Eq. 
(7), we define the second-order derivative of the ground- 
state energy as the anisotropic susceptibility Xa, 

_ d\EjN) _zd%, 

The numerical results for Xa as a function of A arc shown 
in Fig. 2 for the IDC and 2DSL models respectively. It 
is clear from the figure that there is no singular behavior 
for the IDC model. By comparison the exact calcula- 
tion H] gives an essential singularity at A = 1, with the 
result that any finite order of derivative with respect to 
A is indeed continuous. However, the anisotropic sus- 
ceptibility for the 2DSL model clearly shows a singular 
behavior (except for the low-order LSUB2 scheme) . Fur- 
thermore, although the values of the terminating points 
Am are different for the LSUB4 and LSUB6 cases, both 
schemes yield similar critical behavior, namely 

Xa ^ const. X (1 - a;^)"^, x ^ 1, (12) 

with A = 3/2, and x = /S.^/ This clearly suggests that 
Eg/N for the 2DSL model within the LSUBm approxi- 
mation has the following expansion near the terminating 
point, 

^ ^ A„ + B„(l - x2)i/2 + c„(l - x^) 

+Z?„,(1 - x^f'^ + • • • , (13) 

where x = A^/A. 

We note that the spin-wave theory of Anderson |^] 
gives the exponent A = 1/2. It is not difhcult to show 
that our full SUB2 scheme also yields the same value @, 
whereas from the above we see that the LSUBm schemes 
yield A = 3/2. In order to understand this difference, we 
consider the SUB2-m scheme of the 2DSL model with 
m < 14. We observe that the SUB2-m scheme also 
has the singular behavior described by Eq. (13), yielding 
A = 3/2. However, the coefficient Bm ^ as m ^ cxd. 
Therefore the fourth term in Eq. (13) is the dominant 
singular term, yielding A = 1/2 for the full SUB2 scheme 
(i.e., the SUB2-oo scheme). Furthermore, we notice that 
the coefficients Bm also seems to decay by a l/w? rule 
as m increases in the SUB2-m scheme. It seems likely 
that in the LSUBm scheme, the Bm coefficients also de- 
crease as m increases, and the first singular term has the 
coefficient D„i as m — > oo, so that A = 1/2 is recovered. 
Although we only have the LSUB4 and LSUB6 schemes 
for consideration, the results confirm that the value of 



Bm in the LSUB6 scheme is much smaller than that of 
the LSUB4 scheme. Therefore we expect that the exact 
result, namely the LSUBm scheme with m — > oo, has the 
value of A = 1/2. We will present the details of the ex- 
pansion Eq. (13) for our various approximation schemes 
elsewhere. 

While our CCM analysis seems to agree with spin-wave 
theory Q for the critical exponent A, there does seem to 
be a real difference in the corresponding predictions for 
the value of the critical point Ac. We firstly consider 
the SUB2-m scheme. As given previously the termi- 
nating point is Aoo — 0.7984. The terminating point in 
the SUB2-m scheme rapidly approaches this value as m 
increases. Interestingly, the 1/w? rule also seems to fit 
well. We attempt to apply the 1/w? rule for the LSUBm 
scheme also. We thus obtain the predicted critical point 
as Ac w 0.92±0.01. Although this is clearly smaller than 
the value 1 predicted by spin- wave theory , the quoted 
error is merely that of least-squares fit at this level, and 
we cannot yet preclude agreement when higher-order cor- 
rections in the LSUBm scheme (i.e., m > 6) are included 
in the extrapolation. 

Finally, we calculate the staggered (sublattice) magne- 
tization A'F- = {Sj)/s for the XX Z models, where the 
expectation value is taken with respect to the ground 
ket and bra states within the same CCM approximation 
schemes Our results show, as expected, that there is 
clearly a difference between the IDC and 2DSL models. 
For the IDC case, the staggered magnetization shows no 
singular behavior at all, even for the LSUBIO scheme. 
However, for the 2DSL model, the singular behavior is 
clearly seen for all high-order schemes beyond the LSUB2 
scheme. For the 2DSL model, at the isotropic point, 
= 0.8514,0.7648,0.7278 in the LSUB2, LSUB4 and 
LSUB6 schemes respectively. From these numbers we ob- 
tain = 0.68 ±0.01 by extrapolation using a 1/m rule 
(which was found to work well for the SUB2-m scheme) . 
This value is somewhat bigger than the corresponding 
values 0.606 from spin-wave theory |Q, and 0.62 ± 0.02 
from series expansion techniques Q, although it agrees 
well with the best of the corresponding Monte Carlo re- 
sults, which vary between 0.68 ± 0.02 and 0.62 ± 0.04 

In conclusion, our high-order CCM analyses not only 
produce with high accuracy the ground-state properties 
for the spin-i XX Z models, but they also enable us to 
study the the quantum phase transitions in a systematic, 
unbiased manner. They clearly show that the LSUBm 
approximations themselves represent a natural extension 
of perturbation theory. In effect, they comprise a well- 
defined analytical continuation or resummation of the 
perturbation theory results within the context of a nat- 
ural and consistent hierarchy of approximations. More 
significantly, this hierarchy is capable of predicting a pos- 
sible quantum phase transition. Another calculation un- 
der consideration is the energy gap of the excited states 
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near the critical point in the LSUBm scheme. Our result 
for the energy gap in the SUB2 scheme is quite similar 
to that of spin-wave theory, but higher-order multi-spin 
correlations have been proved to be significant . 

Finally, we note that the combination of the COM as a 
theoretical framework and the use of computer-algebraic 
techniques to implement it at high orders of approxima- 
tion, has resulted in a formalism which is capable for 
the 2D spin-lattice models of attaining numerical results 
with a precision comparable to those from the much more 
computationally intensive state-of-the-art Monte Carlo 
simulations. It will be of great interest to apply our 
techniques to similar electron-lattice problems of inter- 
est in high-temperature superconductivity which involve 
vacancies on the lattice (e.g., the Hubbard model), where 
Monte Carlo algorithms are not easily applicable, due to 
the infamous fermion sign problem. 
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FIG. 1. Ground-state energy per spin as a function of A for 
the spin-i XX Z model on the 2DSL. Shown are the numerical 
results of the LSUBm scheme with m = 2, 4, 6 and of a Monte 
Carlo calculation of Ref. [10]. 

FIG. 2. The second-order derivative of the ground-state 
energy per spin with respect to A for the IDC and 2DSL 
models as functions of A. Shown are the results of several 
LSUBm schemes. 
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